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We derive mostly analytically the scaling behavior of the number of nonfrozen and relevant nodes 
in critical Kauffman networks (with two inputs per node) in the thermodynamic limit. By defining 
and analyzing a stochastic process that determines the frozen core we can prove that the mean 
number of nonfrozen nodes scales with the network size N as N^^^ , with only A''^'''^ nonfrozen nodes 
having two nonfrozen inputs. We also show the probability distributions for the numbers of these 
nodes. Using a different stochastic process, we determine the scaling behavior of the number of 
relevant nodes. Their mean number increases for large A'' as N^^'^, and only a finite number of 
relevant nodes have two relevant inputs. It follows that all relevant components apart from a finite 
number are simple loops, and that the mean number and length of attractors increases faster than 
any power law with network size. 

PACS numbers: 89.75.Hc, 05.65.+b, 02.50.Cw 



1. INTRODUCTION 

Random Boolean networks are often used as generic 
models for the dynamics of complex systems of interact- 
ing entities, such as social and economic networks, neural 
networks, and gene or protein interaction networks [1]. 
The simplest and most widely studied of these models 
was introduced in 1969 by Kauffman [2] as a model for 
gene regulation. The system consists of N nodes, each 
of which receives input from K randomly chosen other 
nodes. The network is updated synchronously, the state 
of a node at time step t being a Boolean function of the 
states of the K input nodes at the previous time step, 
t — 1. The Boolean updating functions are randomly as- 
signed to every node in the network, and together with 
the connectivity pattern they define the realization of the 
network. For any initial condition, the network eventu- 
ally settles on a periodic attractor. Of special interest 
are critical networks, which lie at the boundary between 
a frozen phase and a chaotic phase [3, 4]. In the frozen 
phase, a perturbation at one node propagates during one 
time step on an average to less than one node, and the at- 
tractor lengths remain finite in the limit N ^ oo. In the 
chaotic phase, the difference between two almost identi- 
cal states increases exponentially fast, because a pertur- 
bation propagates on an average to more than one node 
during one time step [5]. 

The nodes of a critical network can be classified accord- 
ing to their dynamics on an attractor. First, there are 
nodes that are frozen on the same value on every attrac- 
tor. Such nodes give a constant input to other nodes and 
are otherwise irrelevant. They form the frozen core of the 
network. Second, there are nodes whose outputs go only 
to irrelevant nodes. Though they may fluctuate, they arc 
also classified as irrelevant since they act only as slaves to 
the nodes determining the attractor period. Third, the 
relevant nodes are the nodes whose state is not constant 
and that control at least one relevant node. These nodes 
determine completely the number and period of attrac- 
tors. If only these nodes and the links between them 



are considered, these nodes form loops with possibly ad- 
ditional links and chains of relevant nodes within and 
between loops. The recognition of the relevant elements 
as the only elements influencing the asymptotic dynam- 
ics was an important step in understanding the attrac- 
tors of Kauffman networks. The behavior of the frozen 
core was first studied by Flyvbjerg [6]. Then, in an an- 
alytical study of K = 1 networks Flyvbjerg and Kjaer 
[7] introduced the concept of relevant elements (though 
without using this name). The definition of relevant ele- 
ments that we are using here was given by BastoUa and 
Parisi [8, 9]. They gained insight into the properties of 
the attractors of the critical networks by using numerical 
experiments based on the modular structure of the rele- 
vant elements. Finally, Socolar and Kauffman [10] found 
numerically that for critical K = 2 networks the mean 
number of nonfrozen nodes scales as N^^^, and the mean 
number of relevant nodes scales as N^^^. The same result 
is hidden in the analytical work on attractor numbers by 
Samuelsson and Troein [11], as was shown in [12]. 

In this work, we go a step further by deriving these 
power laws analytically for a more general class of net- 
works, and by showing the asymptotic probability distri- 
bution of nonfrozen and relevant nodes in terms of scal- 
ing variables. We also obtain results for the number of 
nonfrozen nodes with two nonfrozen inputs, and for the 
number of relevant nodes with two relevant inputs. The 
outline of this paper is the following. In the next section 
we define the class of networks that we are investigating. 
In Section 3, we introduce a stochastic process that de- 
termines the frozen core of the network starting from the 
nodes whose outputs are entirely independent of their in- 
puts. Then, in Section 4, we analyze the Langevin and 
Fokker-Planck equations that correspond to this stochas- 
tic process and that lead to the scaling behavior of the 
number of nonfrozen nodes. In order to identify the rel- 
evant nodes among the nonfrozen ones, we introduce in 
Section 5 another stochastic process. This process also 
enables us to find their scaling behavior. Finally, we dis- 
cuss in the last section the implications of our results. 
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TABLE I: The 16 update functions for nodes with 2 inputs. 
The first column lists the 4 possible states of the two inputs, 
the other columns represent one update function each, falling 
into four classes. 



2. CRITICAL K = 2 NETWORKS 

The networks we are studying in this paper are the 
K = 2 critical networks. In these networks each node 
has 2 randomly chosen inputs. The 16 possible update 
functions are shown in table I. 

The update functions fall into four classes [5]. In the 
first class, denoted by J- , are the frozen functions, where 
the output is fixed irrespectively of the input. The class 
Ci contains those functions that depend only on one of 
the two inputs, but not on the other one. The class C2 
contains the remaining canalizing functions, where one 
state of each input fixes the output. The class TZ contains 
the two reversible update functions, where the output is 
changed whenever one of the inputs is changed. Critical 
networks are those where a change in one node propa- 
gates to one other node on an average. A change propa- 
gates with probability 1/2 to a node that has a canalizing 
update function Ci or C2 , with probability zero to a node 
that has a frozen update function, and with probability 
1 to a node that has a reversible update function. Con- 
sequently, if the frozen and reversible update functions 
are chosen with equal probability, the network is critical. 
Usually, only those models are considered where all 16 
update functions receive equal weight. We here consider 
the larger set of models where the frozen and reversible 
update functions are chosen with equal (and nonzero) 
probability, and where the remaining probability is di- 
vided between the Ci and C2 functions. Those networks 
that contain only Ci functions arc different from the re- 
maining ones. Since all nodes respond only to one input, 
the link to the second input can be cut, and we are left 
with a critical K = 1 network, which was already dis- 
cussed in [7, 12, 13] and will not be discussed here. All 
the other models, where the weight of the Ci functions is 
smaller than 1, fall into the same class [12]. The treat- 
ment presented in the following, is based on the existence 
of nodes with frozen functions, and it therefore applies 
to all critical models with a nonzero fraction of frozen 
functions. Networks with only canalyzing functions have 
to be discussed separately. 

Let Nf be the number of nodes with a frozen function, 
Nr the number of nodes with a reversible function and 
A^ci and N^^ the number of nodes with a Ci and a C2 
function. We define the systems we are going to consider 



through parameters a = N^JN, (i = Nr/N = Nf/N, 
7 ~ NC2/N. These parameters give the fraction of each 
type of nodes in the network. In the next two sections, we 
determine the properties of the frozen core in the large N 
limit by starting from the nodes with a frozen function. 



3. A STOCHASTIC PROCESS THAT LEADS TO 
THE FROZEN CORE 

We consider the ensemble of all networks of size N and 
with fixed parameters a, (3, 7. All nodes with a frozen up- 
date function are certainly part of the frozen core. We 
now construct the frozen core by determining stepwise 
all those nodes that become frozen due to the influence 
of a frozen node. In the language of [10], this process 
determines the "clamped" nodes. Initially, we place the 
nodes in four containers labelled !F, Ci, C2, and TZ. These 
containers contain Nf , Nc-^ , Nc^ , and Nr nodes initially. 
Since these numbers change during our stochastic pro- 
cess, we denote the initial values as Nj"-\ iV™% iV™\ 
and A^™*, and the total number of nodes as A^™\ We 
treat the nodes in container Ci as nodes with only one 
input and with the update functions "copy" or "invert" . 
The contents of the containers will change with time. The 
"time" we are defining here is not the real time for the 
dynamics of the system. Instead, it is the time scale for 
a stochastic process that we use to determine the frozen 
core. During one time step, we remove one node from 
the container and determine all those nodes, to which 
this node is an input. A node in container Ci chooses this 
node as an input with probability 1/A^. It then becomes 
a frozen node. We therefore move each node of container 
Ci with probability 1 /N into the container T. A node in 
container C2 chooses the selected frozen node as an input 
with probability 2/N. With probability 1/2, it then be- 
comes frozen, because the frozen node is with probability 
1/2 in the state that fixes the output of a C2-node. If the 
C2-node does not become frozen, it becomes a Ci-node. 
We therefore move each node of container C2 during the 
first time step with probability 1 /N into the container J^, 
and with probability 1/A^ into the container Ci. Finally, 
a node in container TZ chooses the selected frozen node as 
an input with probability 2/N and becomes a Ci-node. 
We therefore move each node of container TZ during the 
first time step with probability 2/N into the container 
Ci. In summary, the total number of nodes, N, decreases 
by one during one time step, since we remove one node 
from container and some nodes move to a different 
container. The removed nodes arc those frozen nodes for 
which we already have determined whose input they arc. 
Then, we take the next frozen node out of container 
and determine its effect on the other nodes. We repeat 
this procedure until we cannot continue because either 
container is empty, or because all the other containers 
are empty. If container J- becomes empty, we are left 
with the nonfrozen nodes. We shall see below that most 
of the remaining nodes are in container Ci , with the pro- 
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portion of nodes left in containers C2 and TZ vanishing in essentially frozen, with only a finite number of nonfrozen 



the limit iV' 



00. Then, the nonfrozen nodes can be nodes in the limit A'^* 



00. If we now choose the in- 



connected to a network by choosing the input (s) to every 
node at random from the other remaining nodes. If all 
containers apart from container T are empty at the end, 
the entire network becomes frozen. This means that the 
dynamics of the network go to the same fixed point for 
all initial conditions. 

Let us first describe this process by deterministic equa- 
tions that neglect fluctuations around the average change 
of the number of nodes in the different containers. As 
long as all containers contain large numbers of nodes, 
these fluctuations are negligible, and the deterministic 
description is appropriate. The average change of the 
node numbers in the containers during one time step is 
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puts for these nodes, we obtain simple loops with trees 
rooted in the loops. This property of the frozen phase 
was also found in [10]. 

For the critical networks that this paper focuses on, we 
have NJ" ~ iV™' = (3N'^™, and the stochastic process 
stops at Nf = 1 = /3N'^/N'''\ This means that 



(3 



(4) 



The number of nonfrozen nodes would scale with the 
square root of the network size if the deterministic ap- 
proximation to the stochastic process was exact. We shall 
see below that including fluctuations changes the expo- 
nent from 1/2 to 2/3. The final number of C2-nodcs for 
the deterministic process for the critical networks is 7//?, 
which is independent of network size, and the final num- 
ber of 7?.-nodes vanishes due to Nr = Nf. We shall see 
below that the fluctuations change these two results to a 
(7V™*)i/3-dependence. 

Introducing n = N/N'"' and Uj = N^/N"'' for j = 
r, /, ci,C2, equations (3) simplify to (using A^,*"* = TV™*) 



The number of nodes in the containers, N, can be used 
instead of the time variable, since it decreases by one dur- 
ing each step. The equation for Nr can then be solved 
by going from a difference equation to a differential equa- 
tion. 



ANr ^ dNr 

AN ~ d/V 
which has the solution 
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Similarly, wc find 
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For Nf"' < Nr"', we obtain Nf = at a nonzero value 
of N, and the number of nonfrozen nodes is proportional 
to N'"'. We are in the chaotic phase. For Nj"' > N^"' , 
the values Nr and Nc2 will sink below 1 when N be- 
comes of the order V N'"' . For smaller N , there are 
only T and Ci nodes left, and the second term contribut- 
ing to Nf and Nc^ in (3) can be neglected compared 
to the first one. When Nf falls below 1, there remain 
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I3n^ = rif 

n — 2/3n^ — •yn'^ 



This means that our stochastic process remains invari- 
ant (in the deterministic approximation) when the initial 
number of nodes in the containers and the time unit are 
all multiplied by the same factor. For small n, the ma- 
jority of nodes are in container Ci, since ric^ ~ n — 0(n^). 
Now, if we choose a sufficiently large N"" , n reaches any 
given small value while Nf = Nr = fin^N'"' is still large 
enough for a deterministic description. We can therefore 
assume that for sufficiently large networks Nf/N = f3n 
becomes small before the effect of the noise becomes im- 
portant. This assumption will simplify our calculations 
below. 



4. THE EFFECT OF FLUCTUATIONS 

The number of nodes in container Ci that choose a 
given frozen node as an input is Poisson distributed with 
a mean Nc^/N and a variance Nc^/N. We now as- 
sume that n is small at the moment where noise be- 
comes important, i.e., that the variance of the noise 
NrjN = Ticjn = I - [2(3 + 7)n = 1 - 0{n) is unity. 
The number of nodes in containers C2 and TZ that choose 
a given frozen node as an input is Poisson distributed 
with a mean and a variance 2{Nc2 + Nr)/N. The fluc- 
tuation around the mean can be neglected as this noise 
term is very small compared to Nr and Nc2, the final 



nodes of type Ci. The network is ^^lues of which are large for sufficiently large N"". We 
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therefore obtain the stochastic version of equations (1) 
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AN,. 



AN = 
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(5) 



The random variable ^ has zero mean and unit variance. 
As long as the nj change little during one time step, we 
can summarize a large number T of time steps into one 
effective time step, with the noise becoming Gaussian 
distributed with zero mean and variance T. Exactly the 
same process would result if we summarized T time steps 
of a process with Gaussian noise of unit variance. For 
this reason, we can choose the random variable ^ to be 
Gaussian distributed with unit variance. 

Compared to the deterministic case, the equations for 
Nr and N^^ are unchanged, and we have again Nr = 
7V2Ar;"7(Ar™)2 and N^^ = Af27v™/(iV™)2. Inserting 
the solution for Nr into the equation for Nf, we obtain 



dN^ 
dN 
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(3N 



N m 
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with the step size dN ~ 1 and (^^) = 1. (In the con- 
tinuum limit dN — > the noise correlation becomes 
{£,{N)£,{N')) = d{N-N')). This is a Langcvin-cquation, 
and we will now derive the corresponding Fokker- Planck- 
equation. Let P{Nf,N) be the probability that there 
are Nf nodes in container TZ at the moment where there 
are N nodes in total in the containers. This probabil- 
ity depends on the initial node number Nini, and on the 
parameter (3. The sum 



We have 

P{Nf,N-l) = 

P{Nf + ANf, N)g{ANf\Nf + ANj, iV)d(A7V/) 



d 



dN 
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(PiNf, N)g{ANf\Nf, N))ANf 
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{P{Nf,N)g{ANf\Nf,N)){ANf)' 
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iP{Nf,N){ANf)) 



2dNj 



iP{Nf,N){iANff)) 



The mean change (ANf) during one step is (ANf) = 

+ -^T^j, and the mean square change is {{ANfY) ~ 1. 
This gives the Fokker-Planck equation for our stochas- 
tic process 

dP d fNf ^ pN\p^ld^ 



dN dNf V N J 
We introduce the variables 



2dNj 



N 



^-^and^- 



(8) 



and the function f{x,y) = {N'"^^ / P{N f , N). We 
will see below that f{x,y) does not depend explicitely 
on the parameters N"^'^ and f3 with this definition. The 
Fokker-Planck equation then becomes 



Let W{N) denote the probability that N nodes are left 
at the moment where Nf reaches the value zero. It is 



poo 

J2 P{Nf,N)c, / P{Nf,N)dNf 

,r__, Jo 



Nf = l 



is the probability that the stochastic process is not yet 
finished, i.e. the probability that Nf has not yet reached 
the value at the moment where the total number of 
nodes in the containers has decreased to the value N. 
Since systems that have reached Nf = are removed 
from the ensemble, we have to impose the absorbing 
boundary condition P{0,N) = 0. Let g{ANf\Nf,N) 
denote the probability that Nf decreases by ANf during 
the next step, given the values of Nf and N. 



W{N) = 
Consequently, 

W{N) = 



P{Nf,N)dNf~ P{Nf,N -l)dNf. 
Jo 



P{Nf,N)dNf 



ON Jo 



d 



(iV™//3)-i/3_ViV fix,y)d2 

(iV™//3)-2/3_y^ f{x,y)dx 
(iV™//3)-2/3G'(y) 



(10) 



with a scaling function G{y). W{N) must be a normal- 
ized function, J^WiN)dN = ]^ G{y)dy = 1. This 
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FIG. 1: The function W(iV)(7V""//3)^/^ vs N/{N'"'/pf^^ for 
P = 0.25 and iV™' = 2^'', 2^°, 2^^ Furthermore, 

the graph contains a curve with /3 = 0.125, A*' = 2^^' and a 
curve with (3 = 0.5, iV = 2". The curves aU collapse, con- 
firming the existence of a scaling function G(y). The dashed 
line is a power law ~ 1/ \/7V. 



condition is independent of the parameters of the model, 
and therefore G{y) and f{x,y) are independent of them, 
too. which justifies om' choice of the prcfactor in the defi- 
nition of f{x, y). By integrating equation (9) over x from 
to infinity and by using /(O, y) = /(oo, y) = we obtain 



0, 



which gives us a second relation between f{x,y) and 
Giy): 



VyG{y) 



1 dl 

2 dx 



(11) 



x=0 



The mean number of nonfrozen nodes is 



N 



NW{N)dN = (iV™//3) 



2/3 



G{y)ydy , 



(12) 

which is proportional to (A^™*//?)^^"^. We did not succeed 
in extracting an explicit expression for the function G{y). 
It can be determined by running the stochastic process 
described by the equations (5) on the computer. The 
result is shown in Figure 1, and an almost perfect fit to 
this result is given by 



Giy) ~ 0.25e-^' /^(l - O.S^y + 3y)/Vv- 



(13) 



For small y, the data show a power law G{y) oc y^^^^. 
We can obtain this power law analytically by solving the 
Fokker-Planck equation (9) in the limit of small y. In 
this limit, the term proportional to y^^^ can be dropped, 
and wc have the simpler equation 



y— + f + - — 
dy 2 dx 



2dx^ 



0. 



(14) 
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FIG. 2: P{Nf,N) vs Nf/VN for N'"' = 2^' and /5 = 1/4 
for different N. The thick solid line is the theoretical re- 
sult Eq. (16), which is approached in the limit of small 



The general solution has the form f{x,y) 
J2u Cuy'^fu{x), with the functions satisfying 



2{,. + l)U + xfl + fJ ^0. 



(15) 



The solution is 



^/2 



C2 iFi 



11x2 

2' 2' y 



with two constants Ci and C2 , and with H denoting the 
Hcrmitian functions, and iFi the appropriate hypergeo- 
mctric functions. We expect / to be analytical in y for 
small y, which means that v = 0,1,2^ ... . For sufficiently 
small y, only the term v ^ contributes, and due to the 
absorbing boundary condition we have C2 = 0. We ob- 
tain therefore for small y 



f{x,y) = coxe 



(16) 



From our numerical result (13), together with (11), we 
find Co = 0.5. Inserting Eq. (16) into Eq. (10), we obtain 
for small 



W{N) 



(3 



-1/3 



Co 



2VN 



(17) 



In Eq. (16), the function f{x,y) is independent of y. 
This means that for sufficiently small N the function 
P{Nf,N) depends only on the ratio Nf/VN. This is 
also confirmed by our computer simulations (see Fig. 2). 

We can obtain a set of solutions of Eq. (9) with the 
Ansatz f{x,y) — '^^y^ fu{z) with z ^ x — y"^^^. The 
resulting equation for f^, is identical to Eq. (15) for 
/i/, which was valid for small y. However, an analyti- 
cal expression for the function G{y) can only be given 
if an expansion of the initial condition P{Nf,N) = 
6{Nf — /37V™*) in terms of known solutions can be found. 
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5. RELEVANT NODES 




iV/('iV""/pj 



FIG. 3: The function W^(ArO(iV'"V/3)^''^ vs Nr/iN'"* /pf^^ 
for (3 = 0.5 and /3 = 0.125 and for iV'"' = 
216 217^ 2i», 2^^, 2^°, 2^^ The 12 curves converge with increas- 
ing N towards an asymptotic curve, confirming the existence 
of an asymptotic scaling function F{s). The dashed line shows 
the function F{s) obtained using the data for G{y) obtained 
from the same simulation and Eq. (18). 



The probability Wr{Nr) that Nr nodes arc left in con- 
tainer TZ at the moment where container J- becomes 
empty, is obtained from the relation 



Defining 



and 



Nr 



(Af™//3)V3 



F{s) 



2yi 



(18) 



and remembering W{N)dN — Wr{Nr)dNr, wc find 

WriNr) = [N'-^'/py^'^F^s) . (19) 
The mean number of nodes left in in container TZ is 



Nr 



Wr{Nr)NrdNr = [N''^'' / (3f ' ^ / sF{s)ds 

Jo 



(7V™V/3)i/3 / y^G{y)dy 



(20) 



The number of nodes left in container C2 is -/Vc2 = 

ii/mr. 

We thus have shown that the number of nonfrozen 
nodes scales with network size N"" as (A^™')2/3^ with 
most of these nodes receiving only one input from other 
nonfrozen nodes. The number of nonfrozen nodes receiv- 
ing two inputs from nonfrozen nodes scales as (A^*"*)i/3. 
We have found scaling functions that describe the proba- 
bility distribution for these two types of nodes in the limit 
of large network size. Our next task will be to connect 
these nonfrozen nodes to a network. This is a reduced 
network, where all frozen nodes have been cut off. 



Let us start from the result obtained from the stochas- 
tic process of the previous two sections. Each time we 
run this process we obtain A'' nonfrozen nodes. Out of 
these, Nj. (A^C2) nodes receive input from two other non- 
frozen nodes and have a reversible (canalizing C2) update 
function. We define the parameter 



Nr + Nc. 



N 



(l+7//3)y'/^ 



(21) 



which has a probability distribution /(a) that is deter- 
mined from the condition f(a)da = G{y)dy, 



/(«) 



G 



3ai/3(i + ^//3)2/3 Ui + -^//3 



2/3N 



(22) 



Just as (?(?/), the function /(a) is the exact probability 
distribution only in the thermodynamic limit A^™' 00. 
We determine the relevant nodes by a stochastic pro- 
cess that removes iteratively nodes that are not relevant. 
Each of the N nonfrozen nodes chooses its input (s) at 
random from the nonfrozen nodes. There are altogether 
iV(l -I- a/y/N) inputs to be chosen, and consequently 
the nonfrozen nodes have together A^(l + a/\/N) out- 
puts. The number of outputs of a node is Poisson dis- 
tributed with the mean value (1 -I- a/^fN). The frac- 
tion exp(— 1 — a/VN) of nodes have no output. They 
are the leaves of the trees of the network of nonfrozen 
nodes, and we therefore know that they are not relevant. 
We put them in container number 1. Their number will 
change during the stochastic process that determines the 
relevant nodes. The other nodes are placed in container 
number 2. Their number is Ni ("labelled"), and it will 
be reduced until only the relevant nodes are left. The 
total number of outputs of the nodes in container 2 is ini- 
tially N{l + a/^/N), while their total number of inputs is 

iV(l-|-a/-\/]V)(l — exp(— 1 — a/-\/]V)). Now, we remove one 
node from container 1 and connect its input (s) at random 
to the outputs of the nodes in container 2. The chosen 
output (s) are cut off. If a node whose output is cut off has 
no other output left, we move the node from container 2 
to container 1. It cannot be a relevant node since rele- 
vant nodes influence other relevant nodes. We iterate this 
procedure, until there is no node left in container 1. The 
nodes remaining in container 2 are the relevant nodes. 
During the entire process, the number of outputs in con- 
tainer 2 is identical to the number of inputs in container 
1 and 2. As long as container 1 is not empty, there are 
more outputs in container 2 than inputs, and only when 
the process is finished do the two numbers become iden- 
tical. We can therefore simplify the stochastic process by 
removing container 1 altogether. We simply have to con- 
tinue cutting of outputs from nodes in container 2 and 
removing nodes with no outputs, until the total number 
of outputs of the nodes in container 2 has become identi- 
cal to their total number of inputs. The remaining nodes 



are relevant, and we have then M^^^"''' = Nrei- These 
nodes can then be connected to a network by connecting 
the inputs and outputs pairwise. 

In order to derive analytical results, it is useful to run 
this process backwards. Starting with N nodes with no 
outputs, adding outputs at random will eventually gener- 
ate the Poisson distribution of the number of outputs per 
node that we have started with. The reverse stochastic 
process is therefore defined by the following rule: Be- 
gin with an empty container (former container 2) and N 
nodes outside the container. Most of these nodes have 
one input, and the fraction a/^/N have two inputs. Add 
an output to a randomly chosen node. Put this node in 
the container. Add another output to a randomly cho- 
sen node (choosing every node with equal probability, 
whether the node is inside or outside the container). If a 
node from outside the container is chosen, put it in the 
container. Eventually, the total number of outputs in 
the container will become larger than the total number 
of inputs in the container. The container contains the 
relevant nodes at the moment when the inputs equal the 
outputs for the last time. 

In order to show that the number of relevant nodes 
scales with a/ZV, we define a scaling variable 



t 



Ni 



During one step, an output is added to nodes that are 
already in the container with probability Ni/N . Let No 
count the number of outputs that have been added to 
nodes in the container, i.e.. No = (total number of outputs 
in the container) ^Ni. Then the average rate of increase 
of No is given for sufficiently large N by 



,dNo. 
'dNi' 



N^ 
N 



Let Ni count the number of nodes in the container with 
two inputs. Their rate of increase is 

dNi a 



dNi 

Consequently, the probability distribution for No is given 
by 



1 



(23) 



and the probability distribution for iV^ is given by 




0,001 



FIG. 4: The function Ci(t) as obtained by running the 
stocliastic process described in this section. The dotted line 
corresponds to the function 0.25t, which is a good fit to Ci{t) 
for small t. 



The stochastic process can be viewed as a random walk 
that steps to the right with a rate t and to the left with 
a rate a. It is finished when Ni = No for the last time, 
i.e. when the walk leaves the origin for the last time. We 
determined the probability distribution Ca (t) for this last 
exit time from the origin by a computer simulation. It 
is shown in Fig. 4 for a = 1. For small t, it increases 
linearly in t, because the probability of making a step to 
the right is proportional to t for small times. For a = 0, 
we can obtain an analytical result from the relation 



Coit) 



dPo{Q,t) 
dt 



^ te 



-t'/2 



(25) 



Since we were able to write the stochastic process in 
terms of t and a alone, the probability distribution for 
the number of relevant nodes depends only on the com- 
bination Nrei/\fN and on the parameter a, 

Pa{Nrel)dNrel = Ca (iVref/Vlv) dNrel/^N- (26) 

The relation between N and a is obtained using Eq. (8) 
and (21): 



N = 



TV*"* 
J+1 



1/3 



Taking into account the probability distribution (22) of 
the parameter a, we obtain the scaling behavior of the 
number of relevant nodes, 

= / ^«^(«)^'' ( (iv4(!3^7))^/0 (^) 

(27) 

The error made by taking the upper limit of the integral 
to infinity vanishes for _/V™* oo. We introduce the 
scaling variable 

AT . 
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z 

FIG. 5: The function P{z) for 7//? = (solid line) and 7//3 = 
4 (dashed line) . The results were obtained by running the two 
coupled stochastic processes for 10^ samples. 

which has then the following probability distribution 

The probability distribution for the number of relevant 
nodes depends for large N^™ only on the scaling vari- 
able z. We determined numerically the function P{z) by 
combining the two stochastic processes described in this 
paper. First, we determined a value of a using the process 
of Section 4. Then, we used this value of a to determine 
the last exit time of the stochastic process of this section, 
giving a value of z. The shape of the curves P{z) depends 
on the value of 7//?, and the results are shown in Fig. 5 
for 7//3 = and 7//? = 4, which is the original Kauffman 
model, where each update function has the same weight. 
It is easy to check analytically that 

lim P(z) = V2^/4(l + 7//3)^/^ . 

z— >0 

The mean number of relevant nodes is 

roc / j^ini \ 1/3 (-oo 

Nrel=J Nrelp{Nrel)dNrel= ij^j J zP{z)dz , 

^ ° (30) 

i.e., it is proportional to (iV™*)-'^/'^. Finally, let us give the 
probability distribution for the number of relevant nodes 
with two relevant inputs. Let m denote the number of 
relevant nodes with two relevant inputs and P(rn; z)dz 
the probability of having the number of relevant nodes 
in the interval [N^gi^z), N^eiiz + dz)], with m of them 
having two relevant inputs. Using Equations (23) and 
(24), we can express P as 

P(™;.) = ^"<i<,{iflc(^) 

Po (m|za-i/3) P, (m|za-i/3) 
''j:iPo{l\za-^/^)P,{l\za-^/^) ■ 

(31) 



As Po and Pi decay exponentially fast with increasing 
m, the mean number of relevant nodes with two inputs 
is finite. 



6. CONCLUSIONS 

In this paper, we have obtained the asymptotic prob- 
ability distributions in the limit of large network size for 
the number of nonfrozen nodes, the number of nonfrozen 
nodes with two nonfrozen inputs, the number of relevant 
nodes, and the number of relevant nodes with two rele- 
vant inputs. The mean values of these quantities scale 
with network size TV™' as a power law in TV™'^ with the 
exponent being 2/3, 1/3, 1/3, and respectively. The 
implications of the results arc manifold. 

First, the notion that these networks are "critical" is 
now corroborated by the existence of power laws and 
scaling functions. Originally, it was expected that the 
quantities that display the scaling behavior should be 
the attractors of the network [2]. In the meantime, it 
has become clear that mean attractor numbers do not 
obey power laws [11]. It is the number of nonfrozen and 
relevant nodes that show scaling behavior. 

Next, let us compare the results to those of critical 
K — 1 networks. A K ~ 1 critical network with N nodes 
corresponds to the nonfrozen part of a critical K ~ 2 
network for a = 0. In this case, the probability distribu- 
tion of the number of relevant nodes is given by Eq. (26) 
with a = 0, 

p,iNr.,).—C„^—j-—e . (32) 

The mean number of relevant nodes is proportional to 
\/N. When these relevant nodes are connected to a net- 
work by pairwise connecting the inputs and outputs, one 
obtains a set of simple loops. From [13], we know that 
there is a mean number of In \/N loops and that the 
number of loops of length / in a critical K = 1 network 
is Poisson distributed with a mean l/l for / <C a/TV- This 
can be easily explained by consindcring the process of 
connecting inputs and outputs: We begin with a given 
node and draw the node that provides its input from all 
possible nodes. Then, we draw the node that provides 
the input to the newly chosen node, etc., until the first 
node is chosen and a loop is formed. For small loop size, 
the probability that the loop is closed after the addition 
of the Zth node is 1/Nrei- Therefore, the probability that 
a given node is on a loop of size I is 1 /Nrei , and the mean 
number of nodes on loops of size Z is 1, and the number 
of loops of length I is Poisson distributed with a mean 
\/l for sufficiently small /. 

Now, the K = 2 critical networks have of the order 
of (iV™' )^/'^ relevant nodes, with only a finite number of 
them having two relevant inputs. The relevant compo- 
nents are constructed from the relevant nodes by pairwise 
connecting inputs and outputs. In the asymptotic limit 
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of very large TV'"* that we are considering, the probabil- 
ity that a randomly chosen relevant node has two inputs 
or two outputs vanishes. Let us again construct a com- 
ponent by starting with one node and choosing its input 
node etc., until the component is finished. If the compo- 
nent is small, it consists almost certainly only of nodes 
with one input and one output and is therefore a simple 
loop. There is no difference between the statistics of the 
small relevant components of a -ftT = 1 critical network, 
and the number of loops of length / is Poisson distributed 
with a mean l/l. The total number of relevant nodes in 
loops of size I < Ic with Ic = e(7V™*)^/^ (with a small e) is 
Ic, and it is a small proportion of all nodes. If there were 
no nodes with two inputs or outputs, the number of com- 
ponents larger than Ic would be [hi Nrei — hilc) = ln(l/e). 
The additional links may reduce this number, which is in 
any case finite. Since these large components contain al- 
most all nodes, they contain almost all relevant nodes 
with two inputs or outputs. 

From these findings, we can obtain results for the at- 
tractors of K = 2 critical networks. The numbers and 
lengths of attractors are determined by the relevant com- 
ponents. We now argue that the mean number and length 
of attractors increases faster than any power law. If we 



remove the components of size larger than Ic and deter- 
mine the mean number and length of attractors for this 
reduced relevant network, we have a lower bound to the 
correct numbers. Now, the reduced relevant network of 
a. K = 2 system is identical to that of a critical X = 1 
system (where the critical loop size is Ic = e\/N)- In 
[13], it was proven that the mean number and length of 
attractors for such a reduced K — 1 system increases 
faster than any power law with network size. We there- 
fore conclude that the same is true for critical K = 2 
networks. 

Earlier, Samuelsson and Troein [11] have derived ana- 
lytically an exact expression for the number of attractors 
of length L of a critical K ~ 2 network in the limit of 
large iV™% and they have pointed out that this implies 
that the mean number of attractors increases faster than 
any power law with _/V™*. Using their calculation, it has 
recently been shown [12] that there is a close relationship 
between K = \ critical networks and the nonfrozen part 
oi K ~ 2 critical networks, and that the results of [11] 
can be most naturally interpreted if the relevant compo- 
nents of these two networks look identical for component 
sizes below the above-given cutoffs. This interpretation 
is placed on a firm foundation by the present paper. 
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